Code
library(edgeR)
library(dplyr)
library(tidyr)
library(ggplot2)
library(scuttle)
library(pheatmap)
library(patchwork)
library(ggbeeswarm)
library(EnhancedVolcano)
library(SingleCellExperiment).volcano <- \(df, title, fdr = 0.05, lfc = 1) {
EnhancedVolcano(df,
x = "logFC", y = "FDR",
FCcutoff = lfc, pCutoff = fdr,
pointSize = 1, raster = TRUE,
title = title, subtitle = NULL,
lab = df[["gene"]], labSize = 2,
drawConnectors = TRUE, widthConnectors = 0.5) +
guides(col = guide_legend(override.aes = list(alpha = 1, size = 3))) +
theme_bw(9) + theme(
aspect.ratio = 1,
legend.title = element_blank(),
panel.grid.minor = element_blank())
}# setup subtype groupings
names(kid) <- kid <- c("Kidney", "RCC")
names(lun) <- lun <- c("Alveoles", "LSCC")
names(tls) <- tls <- c("E_TLS", "SFL_TLS", "PFL_TLS", "Tcell_TLS")
# split by tumor type
dat <- c(list(
RCC = sub[, sub$TumorType == "RCC"],
LSCC = sub[, sub$TumorType == "LSCC"]),
lapply(tls, \(.) sub[, sub$TissueSub == .]))ref <- c(
RCC = unname(kid[1]),
LSCC = unname(lun[1]),
sapply(tls, \(.) "LSCC"))
names(ids) <- ids <- names(dat)
dat <- dat[ids]; ref <- ref[ids]
fit <- lapply(ids, \(.) {
# setup design matrix
df <- data.frame(colData(dat[[.]]))
if (. %in% tls) {
tt <- factor(df$TumorType)
df$TumorType <- relevel(tt, ref[.])
mm <- model.matrix(~0+TumorType, df)
colnames(mm) <- gsub("TumorType", "", colnames(mm))
} else {
st <- droplevels(df$TissueSub)
df$TissueSub <- relevel(st, ref[.])
mm <- model.matrix(~0+TissueSub, df)
colnames(mm) <- gsub("TissueSub", "", colnames(mm))
}
# fit GLM model
dgl <- DGEList(assay(dat[[.]]))
dgl <- calcNormFactors(dgl)
dgl <- estimateDisp(dgl, mm)
fit <- glmQLFit(dgl, mm)
})# setup contrasts
cs <- c(list(
RCC = c(list(TLS = list(kid, tls)), lapply(tls, \(t) list(kid, t))),
LSCC = c(list(TLS = list(lun, tls)), lapply(tls, \(t) list(lun, t)))),
lapply(tls, \(.) setNames(list("RCC"), .)))
cs <- lapply(ids, \(.) {
x <- numeric(ncol(mm <- fit[[.]]$design))
lapply(cs[[.]], \(c) {
if (length(c) == 1) {
i <- match(c, colnames(mm))
x[i] <- 1; x[-i] <- -1
} else {
a <- match(c[[1]], colnames(mm))
b <- match(c[[2]], colnames(mm))
x[a] <- -1/sum(a != 0)
x[b] <- 1/sum(b != 0)
}
return(x)
})
})# run DGE analysis
res <- lapply(ids, \(.) {
lapply(names(cs[[.]]), \(c) {
ht <- glmQLFTest(fit[[.]], contrast = cs[[.]][[c]])
tt <- topTags(ht, n = Inf)$table
data.frame(row.names = NULL,
gene = rownames(tt), tt,
contrast = c, subset = .)
}) |> do.call(what = rbind)
}) |> do.call(what = rbind)
rownames(res) <- NULLtop <- res %>%
filter(!subset %in% tls) %>%
group_by(subset) %>%
filter(FDR < 0.1, contrast == "TLS") %>%
slice_max(logFC, n = 40) %>%
split(.$subset)
for (t in names(top)) {
cat("####", t, "\n")
mtx <- logcounts(sub)[top[[t]]$gene, sub$TumorType == t]
cd <- data.frame(colData(sub))[c("TissueSub")]
hm <- pheatmap(mtx,
main = t, fontsize = 6,
col = rev(hcl.colors(51, "RdBu")),
scale = "row", show_colnames = FALSE, annotation_col = cd)
print(hm); cat("\n\n")
}top <- res %>%
filter(subset %in% tls) %>%
group_by(subset) %>%
slice_max(abs(logFC), n = 40) %>%
split(.$subset)
for (. in names(top)) {
cat("####", ., "\n")
mtx <- logcounts(sub)[top[[.]]$gene, sub$TissueSub == .]
cd <- data.frame(colData(sub))[c("TumorType")]
hm <- pheatmap(mtx,
main = ., fontsize = 6,
col = rev(hcl.colors(51, "RdBu")),
scale = "row", show_colnames = FALSE, annotation_col = cd)
print(hm); cat("\n\n")
}top <- res %>%
filter(!subset %in% tls) %>%
group_by(subset) %>%
filter(FDR < 0.05, contrast == "TLS") %>%
slice_max(logFC, n = 25) %>%
split(.$subset)
for (. in names(top)) {
cat("####", t, "\n")
plt <- ggplot(
filter(gg,
TumorType == .,
gene %in% top[[.]]$gene),
aes(TissueSub, expr, fill = TissueSub)) +
facet_wrap(~ gene, scales = "free_y") +
geom_boxplot(size = 0.1, fill = NA, outlier.color = NA, show.legend = FALSE) +
geom_beeswarm(shape = 21, col = "black", stroke = 0.1, size = 1.2, alpha = 0.8) +
guides(fill = guide_legend(override.aes = list(size = 3, alpha = 1))) +
labs(x = NULL, y = "Expression (logCPM)") +
scale_fill_brewer(palette = "Set2") +
theme_linedraw(9) + theme(
panel.grid = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
strip.background = element_blank(),
legend.key.size = unit(0.5, "lines"),
strip.text = element_text(color = "black", face = "bold"))
print(plt); cat("\n\n")
}top <- res %>%
filter(subset %in% tls) %>%
group_by(subset) %>%
#filter(FDR < 0.05, contrast == "TLS") %>%
slice_max(logFC, n = 25) %>%
split(.$subset)
for (. in names(top)) {
cat("####", ., "\n")
plt <- ggplot(
filter(gg,
TissueSub == .,
gene %in% top[[.]]$gene),
aes(TumorType, expr, fill = TumorType)) +
facet_wrap(~ gene, scales = "free_y") +
geom_boxplot(size = 0.1, fill = NA, outlier.color = NA, show.legend = FALSE) +
geom_beeswarm(shape = 21, col = "black", stroke = 0.1, size = 1.2, alpha = 0.8) +
guides(fill = guide_legend(override.aes = list(size = 3, alpha = 1))) +
labs(x = NULL, y = "Expression (logCPM)") +
scale_fill_brewer(palette = "Set2") +
theme_linedraw(9) + theme(
panel.grid = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
strip.background = element_blank(),
legend.key.size = unit(0.5, "lines"),
strip.text = element_text(color = "black", face = "bold"))
print(plt); cat("\n\n")
}R version 4.3.0 (2023-04-21)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.2.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Europe/Zurich
tzcode source: internal
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] EnhancedVolcano_1.18.0 ggrepel_0.9.3
[3] ggbeeswarm_0.7.2 patchwork_1.1.2
[5] pheatmap_1.0.12 scuttle_1.10.1
[7] SingleCellExperiment_1.22.0 SummarizedExperiment_1.30.1
[9] Biobase_2.60.0 GenomicRanges_1.52.0
[11] GenomeInfoDb_1.36.0 IRanges_2.34.0
[13] S4Vectors_0.38.1 BiocGenerics_0.46.0
[15] MatrixGenerics_1.12.0 matrixStats_0.63.0
[17] ggplot2_3.4.2 tidyr_1.3.0
[19] dplyr_1.1.2 edgeR_3.42.2
[21] limma_3.56.0
loaded via a namespace (and not attached):
[1] beeswarm_0.4.0 gtable_0.3.3
[3] xfun_0.39 htmlwidgets_1.6.2
[5] lattice_0.21-8 Cairo_1.6-0
[7] vctrs_0.6.2 tools_4.3.0
[9] bitops_1.0-7 generics_0.1.3
[11] parallel_4.3.0 tibble_3.2.1
[13] fansi_1.0.4 pkgconfig_2.0.3
[15] Matrix_1.5-4 RColorBrewer_1.1-3
[17] sparseMatrixStats_1.12.0 lifecycle_1.0.3
[19] GenomeInfoDbData_1.2.10 farver_2.1.1
[21] compiler_4.3.0 munsell_0.5.0
[23] codetools_0.2-19 vipor_0.4.5
[25] htmltools_0.5.5 RCurl_1.98-1.12
[27] yaml_2.3.7 pillar_1.9.0
[29] crayon_1.5.2 BiocParallel_1.34.0
[31] DelayedArray_0.26.1 tidyselect_1.2.0
[33] locfit_1.5-9.7 digest_0.6.31
[35] purrr_1.0.1 labeling_0.4.2
[37] splines_4.3.0 fastmap_1.1.1
[39] grid_4.3.0 colorspace_2.1-0
[41] cli_3.6.1 magrittr_2.0.3
[43] S4Arrays_1.0.1 utf8_1.2.3
[45] withr_2.5.0 DelayedMatrixStats_1.22.0
[47] scales_1.2.1 rmarkdown_2.21
[49] XVector_0.40.0 beachmat_2.16.0
[51] evaluate_0.20 ggrastr_1.0.1
[53] knitr_1.42 rlang_1.1.1
[55] Rcpp_1.0.10 glue_1.6.2
[57] rstudioapi_0.14 jsonlite_1.8.4
[59] R6_2.5.1 zlibbioc_1.46.0